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3 I Abstract. The purpose of this paper is to improve the modelization of the rotational mixing which occurs in stellar radiation 

zones, through the combined action of the thermally driven meridional circulation and of the turbulence generated by the shear 

OO ' of differential rotation. The turbulence is assumed to be anisotropic, due to the stratification, with stronger transport in the 

horizontal directions than in the vertical. The main difference with the former treatments by Zahn (1992) and Maeder & Zahn 
(1998) is that we expand here the departures from spherical symmetry to higher order, and include explicitly the differential 
rotation in latitude, to first order. This allows us to treat simultaneously the bulk of a radiation zone and its tachocline(s). 
Moreover, we take fully into account the non-stationarity of the problem, which will enable us to tackle the rapid phases of 
evolution. The system of partial differential equations, which govern the transport of angular momentum, heat and chemical 

_^ I elements, is written in a form which makes it ready to implement in a stellar evolution code. Here the effect of a magnetic field 

VQ . is deliberately ignored; it will be included in forthcoming papers. 

o. 
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Oh 
I 
Q ' In rotating stars, radiation zones are generally not in radiative equilibrium, as was shown by Von Zeipel (1924). Eddington 

(1925) and Vogt (1925) concluded that this thermal imbalance drives a large-scale meridional circulation, which was described 



> 

OO 



1. Introduction - Motivation 



H 



a 



^ in detail by Sweet (1950) for a given rotation law. Mestel (1953) discussed the role of composition inhomogeneities, which build 

. . , up through the nuclear burning of hydrogen; he showed that they would slow down the circulation, and possibly prevent any 

. ^ mixing. Partly for this reason, rotational mixing was largely ignored in further work on stellar evolution. The next significant step 

k> was made by Busse (1982) who showed that angular momentum would be advected in such a way as to achieve a state of zero 

^ circulation, at least in an non-evolving, inviscid star with no stress applied to its surface. 

The subject has been re-examined in Zahn (1992), who assumed that hydrodynamical instabilities due to the shear of differ- 
ential rotation would enforce a rotation law with much larger variation in the vertical than in the horizontal direction. With such 
'shellular' rotation, the problem reduces to one dimension (or rather 1.5 dimension, since departures from spherical symmetry 
are taken in account), and this renders it much more tractable. In particular one can then treat consistently the redistribution of an- 
gular momentum through the meridional circulation, and thus its feedback on the rotation profile, while preserving its advective 
nature. This makes it possible to prove that the boundary condition applied at the surface on the angular momentum flux plays a 
crucial role. In a star which loses angular momentum through a wind, the meridional circulation is enforced by the requirement 
of transporting this momentum to the surface. On the contrary, when the star loses no angular momentum, it relaxes to a state 
where the advection of momentum is balanced by its transport through turbulent motions. Thereafter this treatment was extended 
to inhomogeneous stars, with a general equation of state, and partly adapted to phases of fast stellar evolution (Maeder & Zahn 
1998). 

The first evolutionary sequence calculated according to these prescriptions was that of a 9 M© star (Talon et al. 1997); it was 
followed by extensive grids of models built by Meynet and Maeder (2000), who refined the prescription for the loss of angular 
momentum by taking into account the latitudinal variation of the wind. For massive stars, these models are in better agreement 
with observational data than models without rotational mixing: they account for the observed surface composition, and they 
predict the correct ratio between red and blue supergiants observed in open clusters and in the Magellanic clouds (Maeder & 
Meynet 2001). They also explain the absence of light elements depletion on the blue side of the Li dip, as shown by Talon 
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and Charbonnel (1998). But they fail to correctly reproduce the flat rotation profile in the radiative interior of solar-type stars 
(Mathias & Zahn 1997). In such stars, which are slow rotators, other physical processes must therefore contribute to the transport 
of angular momentum; the most plausible candidates are either magnetic stresses, as advocated already by Mestel (1953), or 
internal gravity waves emitted at the interface with convection zones (Talon et al. 2002; Talon & Charbonnel 2003). 

The purpose of the present paper is to improve the modeling of rotational mixing, in particular during the phases of fast stellar 
evolution, when steep vertical gradients develop, both of angular velocity and of chemical composition. To achieve this, our aim 
is threefold : 

- we extend the linear expansion of all variables to higher order spherical harmonics, in order to capture the tachocline circu- 
lation which is induced in its vicinity by the diff'erential rotation of a convection zone, 

- for the same reason, we treat explicitly the departure from shellular rotation, in the linear approximation, 

- we retain all time derivatives, except those describing the relaxation to hydrostatic and geostrophic balance, in order to better 
resolve the phases of fast stellar evolution. 

In forthcoming papers, we shall introduce a magnetic field, and examine its impact on the thermally driven circulation. 

2. Main assumptions 

Our treatment is based on the assumption that the rotating star is only weakly two-dimensional, and this for two reasons. The first 
is that the rotation rate is sufficiently moderate to allow the centrifugal force to be considered as a perturbation, compared to the 
gravitational field. The second reason rests on the conjecture that the differential rotation induced by the meridional circulation 
gives rise to turbulent motions which are strongly anisotropic due to the stable stratification, with much stronger transport in the 
horizontal directions than in the vertical: i.e. Vi, <*c Vh and Dy <K D/, respectively for the turbulent viscosity and diffusivity. Such 
anisotropic turbulence is observed in the Earth's atmosphere and oceans; in a star we expect it to smooth the horizontal variations 
of angular velocity and of chemical composition, a property we shall invoke to discard certain non-linear terms. The prescriptions 
to be used for these turbulent diff'usivities are discussed and updated in a companion paper (Mathis et al. 2004). 

We thus consider an axisymmetric star, and assume that the horizontal variations of all quantities are small and smooth enough 
to allow their linearization and their expansion in a modest number of spherical harmonics Pi(cos ff). As reference surface, we 
chose either the sphere or the isobar, and write all scalar quantities either as 

X(r, 6) = Xoir) + ^ X,(r)P, (cos 6) (1) 

or 

X(P, 6) = X(P) + ^ Xi(P)Pi (cos 6) . (2) 

l>0 

Let us establish the relation between those two modal expansions. We expand the pressure around the sphere as: 

P(r, 0) = Pair) + ^ P,(r)Pi (cos 6) , (3) 

l>0 

and introduce the radial coordinate of the isobar 

rp(r,9) = r + Y,Ur)P,(cose), (4) 

where r is the mean value of the radius of an isobar Taking the Taylor expansion of P to first order, we have: 

P\r + J] Ur)P,(cos 0), = P^ir) + ^ Pi{r)Pi{cos 0) + I ^ j ^ ^/(r)PKcos 6), (5) 

V Z>0 / />0 \ '' I />o 

and since by definition the pressure is constant on the isobar, we conclude that 

dPo/dr 
If we apply the same procedure to any other variable X, we get 

P;(cos0) (7) 



x\r + Y^^,(r)P,{co^6),e\ = Xo(r) + ^ Xi(r) - (^j^KO 
V />o / ;>() L \ 0/ 



/>o 
and therefore 

X,{r)^X,{r)-{^\pi(r). (8) 

This relation will serve below in §5.1 to calculate the effective gravity. 
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3. Transport of angular momentum 

We start from the momentum equation 

p [d,V + (F ■ V) y] = -VP - pV0 + V ■ ||t|| (9) 

where p is the density, (p the gravitational potential, and ||r|| represents the turbulent stresses. The macroscopic velocity field V is 
the sum of a zonal flow with angular velocity Q.{r, 0) and a meridional flow lAir, 0): 

V = rsm0n(r,e)e^ + tl(r,e). (10) 

The latter can be split into a spherically symmetric part, which represents the contraction or dilation of the star during its 
evolution, plus the meridional circulation 

^l ^re,. + ^^M(r,e), (ii) 

which we expand in spherical functions 



l>0 



dP,(cos 0) 
Ui(r)P, (cos 0) er + V,(r) '\ ' eg 

d0 



(12) 



Using the continuity equation in the anelastic approximation, i.e. V ■ (p1Jm) - 0, we obtain the following relation between t//(r) 
and V,(r): 

Making again use of the continuity equation, the azimuthal component of (|5) can be written as an advection/diffusion equation 
for the angular momentum: 

d sin^ 1 

p— (r^ sin^ 0^) + V ■ (pr^ sin^ 0Q^(M) = — ^5,. (pv„r'*5,o) + dg (pvi, sin^ deO) , (14) 

(jf V / V J fZ \ '' sin ^ '' 

where as in Zahn (1992) we assume that the eflfect of the turbulent stresses on the large scale flow are adequately described by an 
anisotropic eddy viscosity, whose components are v,, and v/, respectively in the vertical and horizontal directions. As discussed in 
Mathis et al. (2004), they act to reduce the cause of turbulence, as observed in laboratory experiments of rotating flows, namely 
here the vertical and horizontal gradients of angular velocity. This explains why the respective fluxes of angular momentum 
contain only these gradients, in contrast with the treatment of Kippenhahn (1963), who considered the effect of an imposed 
anisotropic turbulence, due to thermal convection; in that case the rotation becomes non-uniform, as is actually observed in the 
solar convection zone. Note the introduction of the lagrangian time-derivative d/dt, meaning that the radial coordinate r is the 
mean radius of the layer (either sphere or isobar) which encloses the mass M, , with dMr - Anpp-dr. 
The form of eq. MAI incites to expand the angular velocity as 

Q(r, 0) = Q(r) + Q.ir, 0) = Q(r) + ^ fi;(r) Qi{0), (15) 

with the horizontal average being defined as 

_ C Q.{r,0)s.\rv'0d0 

j^ sin^ede 

and where the horizontal functions Qi{0) satisfy the orthogonality condition 

CQ,{0)sm^0d0 

° ,. , =0. (17) 

/q sin^ede 

These horizontal functions are readily identified: 

r^/(cos0)sin^6id6l i 

e,(0) = f/(cos0)-/, with I,^'-^ — - — =5/,o--5,,2; (18) 

]^ sin^^ede 5 

except for / = 2, these functions Qi reduce to the Legendre polynomials. 
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3.1. Vertical transport of angular momentum 

Taking the horizontal average of equation J14> . and using the assumption that Q. » Q/, we obtain the following vertical advec- 
tion/difFusion equation for the mean angular velocity Q.: 

P^y^) = ^dr (pr^Qf/z) + ^d,. {pvy^,^ . (19) 

Note that only the I - 2 component of the circulation is able to advect a net amount of angular momentum; the higher order 
components of 11m (for instance those induced in its tachocline by a differentially rotating convection zone) do not contribute to 
the vertical transport of angular momentum, as was explained in Spiegel and Zahn (1992). 

3.2. Horizontal transport of angular momentum 

We establish the equation governing the horizontal transport of angular momentum by multiplying eq. ( I19> through sin^ and 
subtracting it from the original form (I14t : 

p— (r^ sin^ 6»q) + V ■ (pr^ sin^ 0Q.%Im) + ^^5, (d/Q U2] = ^^5, (pv,/5,Q) + dg (pvh sm^edgh] ; (20) 

at ^ / v / 5^2 \ J fZ \ ' smO ^ ' 



in the advection term we have again neglected the fluctuation Q. compared to the mean Q. 

The next step is to replace Q(r, B) by its expansion ( I15t in the horizontal functions Q\(8). For I - 2, the equation separates 
neatly into 

p- (r^Qz) - 2plir \lY2(r) - a(r)f/2] = -^5, (pv^r^S^Qj) - 10pv,,Q2, (21) 

Id/ X idlnfr^n) 

with Y^^-—-[p?U^ and a = t— 77 -, (22) 

opr dr ^ ' 2 d In r 

which we can simplify by assuming that the turbulent transport is much more efficient in the horizontal than in the vertical 
direction, i.e. v,, <K v/,: 

p^ (r^Qz) - Ip^r {1Y2 - af/2] = -lOpv/.Qz. (23) 

at ^ ' 

In the asymptotic regime t » r^/v/,, we retrieve the equation for Q2 given by Zahn (1992): 

Vh^iir) = ^r [IViir) - a(r)U2(r)] Q.{r). (24) 

For / > 2 the situation is more intricate, with couplings between terms of dififerent /. However, when making the three 
following assumptions: (i) stationarity {t » r^ Ivh), (ii) strong anisotropy of the turbulent diffusion (v,, « v^), (iii) vertical 
advection negligible compared to the horizontal one (thin layer approximation), the /-components separate, and one recovers the 
result given by Spiegel and Zahn (1992) in their treatment of the solar tachocUne: 

v/,Q,(r) = Qry,(r), (25) 

namely that horizontal diffusion balances horizontal advection. 

4. Transport of chemicals 

The transport of chemicals species is governed by an advection/diffusion equation similar to that for the transport of angular 
momentum: 

pd,Ci +p%l-Vci^ —dr (r^pD„drCi) + . . de {pDh sin 6 dgCj) , (26) 

r^ ^ ' r-^ sm 9 

where c, is the concentration of a given element, and where we assume again that the turbulent difiiisivity is anisotropic. For 
simplicity, we have not included here the term representing the production or destruction of the species due to nuclear reactions, 
which can be easily added. We follow the same procedure as for the angular momentum: we split the concentration in its 
horizontal average and its fluctuation 

Ci{r,0)^ci(r) + c',(r,0), (27) 
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and explicit the velocity field in 

^l^err + ^lM+ BrUf^ (28) 

where, to the radial displacement accompanying the evolution of the star and to the meridional circulation, we add the microscopic 
diffusion velocity V^^^ , which depends on the considered species. We insert these expressions into ( I26> . which leads us to 

p- (cj + c\) + p{nJM.r + Uf)d,-i +p%Im-'Vc'.^ -dr (r^pD,dr (cj + c'.)) + -^--dg (pDh sin edgc'^ . (29) 

at H ^ ' r^smO 

Then we take its horizontal average, making use of the anelastic continuity equation V ■ (p 11m) - 0; we thus obtain the equation 
governing the evolution of the mean concentration cj: 

p^ci + \dr \r^p < c] %iM,r >] + ^5, ^ p-^iUf] = ^3, ir^pD.drl]) , with < / >= ^-^^:^ . (30) 

df r^ L i yi I i yi \ I j^ sinflde 

The subtraction of J30> from ( I29> yields the following equation for the fluctuations of the concentration: 

p—c'-^ + p1lMjdrCl + p%Im ■'^c'. - —dr\r^p<c\1lM,r >] = -^drir^ pDydrc') + , . dg (pDh sin edgc'i). (31) 

Expanding c'. and I/m,/- in Legendre polynomials: 

c'i^Yj^ij(r)Pi(cose) and IIm.,- ^ J]Ui(r)P,(cos9), (32) 

and assuming that the gradient of the concentration fluctuations is negligible compared to that of the mean concentration, i.e. 
|VcJ| <K dr~i, this equation can be recast into 

P ■T-c'iJ + pUldrCj = —dr {r^ pD ,.8 r'Ci J) 2 PDhCi,!- (33) 

As in Chaboyer and Zahn (1992), we further assume that the horizontal diffusivity is much stronger than the vertical one, 
more precisely that Di, » D^,{lhllvf', where //,(/,) is the distance over which c'. changes significantly in the horizontal (vertical) 
direction, and we finally obtain: 

d- ,_ /(/+1) _ 

—Cij + Uidi-c;^ — DhCij. (34) 

at r^ 

For t » r^ jDjj there is an asymptotic state where 

^•' = -/(7TT)Dl^'^^^ ^^^^ 



and then 



(36) 



However, in order to better resolve the fast evolution phases, one may prefer to retain the time derivative in J34I I. and solve the 
equation for the mean concentration J30> in the form: 



d_ 1 



p — Ci + —dr 



2 V '-'.' 



rpl^ 



CijUi 



,>o(2^+l) 



+ \dr [rV^f/f'*] = \dr [r^pD,.dr^] . (37) 



Finally, let us recall how the molecular weight p. is related to the concentrations c,: 

p Y ^'' 

where A, is the mass number (protons+neutrons) and Z, the number of electrons. From ( I34> we can thus derive the following 
advection/diffusion equation for the horizontal fluctuation of the molecular weight A/ = pi/Jl'. 

d , dlniJ . U, 1(1+1) , 

at at Hp r^ 

where we have introduced the pressure scale-height Hp - |dr/dlnP| and the logarithmic gradient V^ = dln/7/dlnP. This 
equation will play a key role in the derivation of the meridional circulation, which we shall discuss in §6. 
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5. Structural properties of the differentially rotating star 

5.1. Baroclinic relation 

In our case where we take a non-uniform and a non-cylindrical rotation law, the centrifugal force does not derive from a potential. 
We are in the baroclinic configuration, where the isobars and the surfaces of constant density do not coincide anymore. To find 
how the density varies on an isobar, we start from the hydrostatic equation: 



- VP = g = -V0 + fc with fc = :r^^V (r^ sin^ o) 



(40) 



where ip is the gravitational potential, and where the local effective gravity g includes the centrifugal force Tc- Taking the curl 
of this equation, we get 



1 11, , 

- -5-Vp A VP = — Vp A g = -V(Q.Y A V(r sin6ir, 
P P 2 

which to first order takes the form 

P_^ ^ [5r(ilV cos 61 sin 6»-5fl(Q^)sin^ el e^, 

P ^ J 

with p' being the variation of the density on the isobar, which we write as 
p(r,0) = 2p,(r)P,(cos0). 

Likewise we expand Q^ in spherical harmonics; since (cf. llSlFTSl 
£l(r, 0) = Q(r) + Y^ Q,(r) (P,(cos 0) - h) 

we have 



;>o 



a\r,0)^ 



Q -2^^2/2 



+ 2Q^Q,P/(cos6i) 



(41) 



(42) 



(43) 



(44) 



(45) 



where we kept only the terms linear in Q/. We expect the largest departures from shellular rotation to occur in the tachoclines, 
where they are enforced by the differential rotation of the adjacent convection zone; judging by the solar case, such departures 
should be rather small, of the order of 1/10, thus justifying the linear approximation. 

Inserting these expansions in (I42> . we reach after some algebra the following expression for the modal amplitudes of the 
relative density fluctuation on an isobar 



^ - '^mr). 



where g is the horizontal average, on the isobar, of the modulus of g, and 



(46) 



^2, 



DM = N?\rdr\n{r)- 2Q.(r)Q.2{r)l 



1 



2>Nl 



■q6i2 



+ 2r 2 5. (n(r)Q,(r)) -^ fA° (-C°_i(5,,,_2 + O^-it^/,,) + B° (-C°^,5,,, + Z)?+i<J/..+2)] 

s>0 ^^ ) 



(47) 



All numerical coefficients involved (Aj', B^, Cj*, Dj', Gj', //", Nf) are given in Appendix A. We shall display here the result for 
I = 2,4, where we keep only the first term Q2 of the expansion of Q: 



^2 = T 

3 



rdr^ 



— [r5,(QQ2)] + -Q02 



D4 = — ^rdr [0.^2) - 2002] . 



(48) 
(49) 



For Q2 = we recover the expression given in Zahn (1992): p2/p - (r /3'g)dr^ . 

This baroclinic equation ( 1461 - 147> plays a key role in linking the density fluctuation on an isobar with the rotation profile. 
It will allow us to close the system formed by the equation for the transport of angular momentum, that for the transport of the 
chemical species and that for the transport of heat, which we shall establish in §6.. 
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5.2. Effective gravity 

Next we examine the redistribution of masses by the centrifugal force and its effect on the gravity inside the star As all other 
scalars, we expand the gravity as: 



g(P, 9) = g{P) + ^ gi{P)Pi{cos 0). 



(50) 



/>o 



The goal of this section is to determine the amplitude of the fluctuation on an isobar, 'gi. The first step will be to calculate the 
perturbation of the gravitational potential (p on the sphere of radius r, and thus to calculate the functions (piir), where we have 
expanded (p as: 



<P(r, 0) = ^^^\r) + (P^'Hr, 6) = Ur) + Yj 4>i{r)Pi{co& 9) with Ur) 



GM(r) 



(51) 



;>o 



We follow the method of linearization developed in Sweet (1950), and expand the pressure and the density around the sphere in 
the same way as (p: 



P(r, 0) = P^^\r) + P^^\r, 0) = Po{r) + ^ Pi(r)P,(cos 0) 

l>0 

p(r, 8) = p^''\r) + p^'Hr, 8) = po(r) + ^pKr)P,(cos 0). 



l>0 

Then, we take the hydrostatic equation (using the classical definition of a potential, contrary to Zahn 1992): 

Vf 1 

= -VS + Tc where Tc = -O^V(r^ sin^ 0), 

P 2 

which we expand to first order: 

Vp(0) ^ _p(O)v0(O) 

vp(i) = -pmv0(i> -p(i'v,^(°' +p'"^rc- 

We eliminate the pressure fluctuation by taking the curl of the latter equation, which gives us: 

^0 go^ ^ '^ 

Next we insert the modal expansion of p*'' and those of the components of the centrifugal force: 

TcAr, e)^Yj «'('')^'('^os 0) Tc,6{r, 0) = - ^ /7,(r)5eP,(cos 0); 

after integration in 0, this yields the modal amplitude of the density fluctuation over the sphere: 



Piir) = — 
§0 



— — 0i +pofl/ + — (rpcfii) 
ar dr 



It remains to insert this expression in the perturbed Poisson equation V^0; - 4nGpi to retrieve Sweet's result: 






1(1 + I)- 47rG dpo- 47tG 



01 ■ 



0/ 



The functions fl;(r) and bi(r) are given by 



poai + — (rpobi) 
ar 



a,(r) = -r[QV)-2Q(r)n2(r)/2 



{^1,0 - 812) 



+ 2rQ.{r) Y, ^s(r) [-^ {C° {olXhs ' hI,NI,6,,,^2) - D^ {clX.ih.^^i - H^Xh^^} 
j>o L^/Vj 



bi(r) = ir[QV)-2Q(r)n2(r)/2 



Sl,2 






(52) 
(53) 

(54) 

(55) 
(56) 

(57) 

(58) 

(59) 

(60) 
(61) 

(62) 
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where we have used the expansion ( 145 > of Q^. The numerical coefficients are given in Appendix A; here we shall quote only the 
results for 1-2,4, with the expansion of Q limited to the term Q2: 



2 -2 

flo = -rQ. bo - 

2 -2 24 - 1 _2 8 - 

fl2 — ——rQ. + ^^rQ£l2 b2 — —rQ. + ■:^rQQ2 



3 

24 - 
04 - -—rQQ2 



35 



35 



b4 = :rTrQQ2 



(63) 
(64) 
(65) 



With this multipolar expansion of the gravitational potential, we are ready to calculate the horizontal fluctuation of the 
effective gravity: 



g^-V(f, + Tc^ -V0o(r) - Yj [V<^/(r)P,(cos 6) + <p,(r)V (P,(cos 6))] + Tc- 



(66) 



Taking the scalar product of g by itself, and retaining only the first order terms, we get the following expansion for g^ over 
the sphere of radius r. 



g\r,e) 



Ar 
Comparing with 



#()('•) \ ^A<t>i){r) 

2 — : flo(r) 



Ar 



+ 2g, 



^(AUr) 



• aiir) 



P/(cos 6*). 



g\r,e) 



gi)' 



ir) + ^?,(r)PKcos 0) = gl(r) + Igair) Y^-giir)?! (cos 0) 



li=0 



1*0 



we identify: 

gi(r) = —: ai(r). 

dr 

Making use of (|8j we get likewise for the variation along an isobar: 



d?- go [po) 



+ ai. 

dr 



It remains to replace Pi by its expression drawn from the 0-component of the hydrostatic equation 

(pi — rbi 
to finally obtain the variation of the effective gravity along an isobar: 



£1 

PC) 



— :jrbi + —ai 

dr gi go 



dr[go, ' 



Keeping only the Q2 term and carrying the expansion to / = 4, we have: 



g dr^2\3 35 / go\3 35 / dr 



(67) 



(68) 



(69) 



(70) 



(71) 



(72) 



(73) 



ll.JJl^(l.TiQ^,J-('±Qn^.±ih 

g dr gl\35 '1 go\35 'j dr[go) 



(74) 



For strictly sheUular rotation (Q2 — 0) we retrieve the expression of Zahn (1992) (after changing the sign of (f>2)' 

^2^ 



g 3 dr\goj dr[gQ^ 



(75) 



The two terms represent respectively the contribution of the centrifugal force and that of the modified mass distribution. The 
latter is obtained by integrating the perturbed Poisson equation ( I60> with the appropriate boundary conditions (see §7). 
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6. Thermal imbalance and transport of heat 

Following closely the procedure described in Zahn (1992), we now undertake to link the meridional circulation with the rotation 
profile and the inhomogeneities of chemical composition. We start from the equation stating the conservation of thermal energy: 



pT 



f.«.vs 



= V-(rVr)+pe-V-f,, 



(76) 



where S is the entropy per unit mass, x the thermal conductivity and e the nuclear energy production rate per unit mass. As 
Maeder and Zahn (1998), we include the flux f /, carried by horizontal turbulence. 

We project this equation on the base of spherical harmonics, starting with the left hand side: 



pT 



dt 



__d5 _- v^ 

/>o 



AS I ^^ dS 
dt dr 



Pi (cos 6) , 



(77) 



where we have used again the expression of the velocity field ( I28> "U - 'e'rr + 1(m, together with the modal expansion of the 
meridional flow (I12> . Recognizing in the horizontal average the production of energy related with the contraction or dilatation of 
the star during its evolution, i.e. TdS /dt - -egra,,, we rewrite ( I76> as 



??z 



/>{) L 



dSi ^, dS 

+ Ui — 

df dr 



P, (cos 0) = V ■ (xVT) +pe + -pig,a, -V-Fh, 



(78) 



thus introducing the lagrangian time derivative. Note that by construction the right hand side has now a zero horizontal average. 
Then we expand the temperature in 



T(P,0) = T{P) + 2 T,(P)Pi(cosd), 



/>o 



from which we deduce the gradient 



VT ^p 



dT 



dT, 



VP 



Ui \— ' Ui / v/^ \— I ~_ 



p 



(79) 



(80) 



/>o 



Multiplying by the conductivity ;t' and taking the divergence, we obtain 



V ■ Crvr) = px 



dT 
dP 



v^ dr, / 



P 



+ V p^ 



dr 

dP 



/>() 



Zai I 
-^P,(cos0) 



VP 

P 



/>() />0 

Next we replace 



^g^-V(f, + rc with \g\ = g{P) + y g,iP)P,icos 0) and Tc = t^^^ V(r sin 0f 

p ^—^ 2 

^ />o 



P 
and expand 



PX^Tx + Yj P^iPiicos 0), V^^ = 47tGp + AnG ^ p/P,(cos 0), V-Tc^fc + Yj /c,/-P/(cos 0). 



(81) 



(82) 



(83) 



/>() 



/>() 



/>o 



The functions /q and fc,i are easily derived. Expressing the components of the centrifugal force in terms of ai and bi introduced 
in ( I58> . we get 



^■^c=i:Uc^.(.%w).^^, 



(r) 



1=0 

and therefore: 



P/(C0S 0) 



fc = —^r ('■^«o) and fcj = —dr (r^fl/) + /(/ + 1)— . 



(84) 



(85) 



Then we explicit the turbulent flux; with the assumption that the turbulent diffusivity is much stronger in the horizontal than 
in the vertical direction (D/, » D,,), its divergence reduces to 



V ■ f ,, = - 



1 
r^ sin 



de[sm0DhpTdgSir,0)] 



(86) 
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which we expand in 



1(1+ I) — ~ 



/>() 



After some arrangements, and keeping only first order terms, this leads us to 



_ d (_AT 



(XVT) +p£ + p6g-V- Fh =< p^— (-4;rGp + fc) + p-^ i 

Z 



dP 



/+p(e+%)> 



(87) 



(88) 






"'^dpK'dp-'^'^dP 



dP r dP 



1 rr-T r\ t /-\ T' 

[-AnGp + fc) + PX-^ [-AnGpi + fcj) + p— | PX' 
_2 /(/ + l) _y , - d f_dT) 



2^^/ 



?2 +pe, - l^l±llpTDi,S,\P,(cos 6). 



As mentioned above, the horizontal average < ... > is zero on a level surface, which means that the star is in average radiative 
equilibrium, including the production of nuclear and gravitational energy. 

It remains to replace all tilded quantities in terms of two only (for each /-component). Here we choose these two variables to 
be the horizontal variations of temperature and molecular weight, instead of Zahn (1992) and Maeder & Zahn (1998) who took 
the density variation. We define 



T ^i p 

and use the equation of state in differential form 
dp _ dP dr dp 

— -"T'^T^'^J 

with the straightforward definitions 6 — —(d Inp/cJ In T)p^^ etc., to express 

0, = ipK, - 5T,. 

We operate likewise with the other variables: 

- = Xt'Vi +X^^i, - = ^T'i'i + e;.A/, 

X e 



and introduce 
K ^ 



M(P) 



L{P) 



and /f 



pC/ -" inr^' - M(P) -■- - (^ + ^^,.„,| 

Switching to r as independent variable, with dP = -p^dr, we cast ( I88> in its final form: 
V • (x'VT) +pe + 'pigrav -^ -Fh 



Y^-T,P,(cos6) 

l>0 



with 



1 - 



fc i^+^gra.) 



AnGp 



g A-nGp A-nGp 



- ~dr{-HTdr'Vl + (1 - 6+Xt)'¥, + (^ +X,)^) - ^ {-HTdr'¥l + {l-6 +XT)'i'l + (f +Xf')^l) 



(e + eg™,.) , . {e + egn.) p„, /(/ + \)Ht 
((/f er - fe6)'V, + {f,£^ + f,ip)A,j - {-6'¥, + tpA,) - "^ 



p 3r 



"¥,+ 



\K ICp 



Here we have used the property that 



PX 



L(r) 



dP g^^^g-dS jg^^^[-V^cf> + V-rc]dV 47TGM(r)-r2ao{ry 



(89) 



(90) 



(91) 



(92) 



(93) 



(94) 



(95) 



(96) 
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To lowest order r^oo - (2/3)r^Q. , according to J61> . which justifies the approximation 
Lir) 



—dT ^v , 

The last step is to put in final form the left hand side of the heat equation j78> . In a medium of varying composition the entropy 
of mixing must be taken into account (Maeder & Zahn 1998). In the simplest case, where the stellar material can be approximated 
by a simple mixture of hydrogen and helium with a fixed abundance of metals, the entropy of mixing can be expressed in terms 
of the molecular weight only; then we have 



dS = Cn 



dT dP du 

T P jj. 



(98) 



where V^d is the adiabatic gradient and O is a function of the metal mass fraction and of ^, the mean molecular weight (see 
Maeder & Zahn 1998). Applying this relation to the fluctuation of entropy around an isobar leads to: 



In the same way, we get for the mean entropy gradient: 
C 



drS = -f (V,d - V - (BV^) 

tin 



(99) 



(100) 



with: V = d In r/d In P and V^ = d In ;u/d In P. 

We refer back to ( I78t to express the left hand side of the thermal energy equation, and take into account that A/ obeys a 
similar advection/diffusion equation ( I39> : this permits the elimination of dA;/df on the l.h.s. and leads us to 



TC„ 



dT, ^dln;;!^ 
dt dt 



Uijr) 



(^ad - V) 



^^'(^>' 



where we cast Ti(r) ( I95l l in its final form, having replaced S i by ( I99l l on the right hand side: 



r, = 2 



1 - 



47tGp 



g AnGp AnGp 



l; + 



(101) 



(102) 



'-d,-{Hrd,^, -{1-5 +xt)^, - (^ +;r.)A/) - ^^LlIlHl (i + |i)vi,, 



I ^ "•" ^gravj 



{{Hrd,^, -(1-6 +XT)'i'i - (if +Xt.)^i) + iUr - fe6 + 6)^, + iU, + f,ip - ip)K,^ 



We recall that /^ + 2 fcjPiicos 0) is the divergence of the centrifugal force (cf. eq.|82]); the functions fQ and fcj are given in 

This expression is similar to that derived in Maeder & Zahn (1998), where ©2 = 'pil'p was used as dependent variable instead 
of *?; - Ti/T here. However we note an important difference, namely that the l.h.s. involves here the only the thermal part of 
the subadiabatic gradient (Vad - V), whereas in Maeder & Zahn (1998) it was (Vo^ - V + ((/)/5)V^j. Also, here the highest order 
derivative at the r.h.s. operates only on ^'i, whereas in Maeder & Zahn (1998) it applied both on ©2 and on A2. We expect 
therefore the present form to be less sensitive, numerically, to steep composition gradients. 

The temperature fluctuation on an isobar is thus governed by an advection/diffusion equation, from which we can derive the 
radial component of the meridional circulation. It allows us to link the circulation with its causes, namely both the rotation profile 
and that of chemical composition. The fact that this equation has been derived for a general law of rotation allows us to treat 
simultaneously the bulk of the radiation zones and the tachoclines. 

The system of equations is now completed: we have an advection/diffusion equation for the transport of the angular momen- 
tum (mean and fluctuating), another for the transport of chemical elements (mean and fluctuating), another for the temperature 
(mean and fluctuating), and we have the baroclinic relation which allows us to close the system. It remains to specify the bound- 
aries conditions of this system, to be applied on the limits of the radiation zone. 



7. Boundary conditions 

We shall now review the differential equations we have established for the various quantities which enter in the problem, and 
state for each of them the appropriate boundary conditions, for a given radiation zone. To be specific, we consider a star with a 
radiation zone located between a convective core and an upper convection zone, and designate by r/, and r, the radius respectively 
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of the base and of the top of that radiative zone. Of course, in a solar-type main-sequence star we have r^ = 0, whereas for a 
massive main sequence star r, - R,R being the radius of the star. 

The equation for the mean angular velocity Q ( I19> is of second order in r, and therefore it requires two boundary conditions. 
They are obtained by evaluating the budget of angular momentum in each adjacent convective zone: 



d 
At 



f 

Jo 



r'^pQ.dr 



1 4 - d 

= -r pQ.U2 atr^Kb — 






r^pQdr 



- --r'^pQ.Uj + Ta at r = r,. 



(103) 



where T^ is the angular momentum flux which is lost at the surface by the stellar wind. (The perturbation Q2 obeys an evolution 
equation which does not include any derivative in r, at least with the approximation we made; therefore it needs no boundary 
condition.) 

Likewise, the equation for the mean concentration J30> is of second order in r, and its two boundary conditions are obtained 
by calculating the evolution in time of that mean concentration c, in the two adjacent convection zones (cf. Palacios et al. 2001). 
We thus have respectively for r, and r;,: 



d 


' 





r, 


dr 




d 







<'i 


At 





f 

J r, 

r 

Jo 



r pdr 



r^pdr 



rV (t/f **'c,) - rV (D„ + Deff ) 5,-c,- - m at r = r„ 
. (f/f «'c,) + rV (D, + Deff) d,ci at r = r^. 



-rVl 



(104) 



where M is the rate of mass loss through the wind. 

Let us examine now the heat equation (I101ll02t . which is of second order in r for the variable *P/. The variable A/ need not to 
be considered here, since it is determined through its evolution equation ( I39> . The boundary conditions on *}*/ are deduced from 
the baroclinic equation J46I47I I. which involves Q and Q2, and their first order derivatives. We must therefore specify all these 
functions at the boundaries of the radiation zone. For Q this is akeady done (ea. ll03t . and for Q/(r) we requke the fluctuations 
to be continuous at the boundaries: 



Q;(r) — Cli^h at r = r^,, i^i(r) — Q/,/ at r = r,, 



(105) 



thus assuming that there are no stresses between the radiative and the convective zone. From the baroclinic relation (I42> we 
deduce that the gradient of the angular velocity must also be continuous, and therefore that 



dr^(r) — dr^b at r = r^, 5fQ(r) = 5,X2, at r = r,, 

dr0.i{r) — dr0.ib at r = r^, dr^i(r) — dr^ij at r = r,. 



(106) 



Since we do not solve here for the rotation law in the convective regions (a formidable task!), the best we can afford is to 
apply a horizontal profile which is inspired from helioseismology, for the base of the outer convection zone. Since the rotation 
rate Q.{r, 0) seems to be independent of r in the solar convection zone, we can take 



5,Q(r) = at r = r,, 5,Q,(r) = at r ^ r,. 



(107) 



For the condition at the boundary of a convective core, we will have to rely on numerical simulations (see the recent results of 
Browning et al. 2004). In the case where we stop the expansion of rotation at the Q.2 term, and where we assume that the vertical 
gradient of Q. is zero, we obtain: 



^A2-5*P2 = -=^^2, 
^ g 

12 r- 
^A4-(5>P4 = - — -nD.2. 
35 g 



(108) 



Finally, let us state the boundary conditions for the perturbed Poisson equation. Unlike the precedent equations above, this 
second order equation must be integrated over the whole star We require regularity at origin, and continuity with a potential field 
at the surface {r - R): 



*, = Oatr = 0, (I + l)d)i + —(pi ^ at r ^ R. 

dr 



(109) 
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8. Conclusion 

The work presented here marks another step in the description of rotational mixing in stellar radiation zones, after the papers by 
Zahn (1992) and Maeder & Zahn (1998). When applied to massive stars, which are fast rotators, such mixing yields much better 
agreement between models and observations, as was discussed by Meynet and Maeder (2000). However numerical problems 
arise during the advanced stages of evolution, when steep gradients of composition and rotation develop. For this reason, our 
main purpose here was to increase the accuracy of the modelization in latitude by adding higher order terms to the expansion in 
spherical harmonics, and also to better resolve the phases of rapid evolution by retaining all time derivatives, except those related 
to the dynamical relaxation. 

This allows us to include the transport occurring in the tachoclines, where the differential rotation imposed by the adjacent 
convection zone generates an octupolar circulation, with no net advection of angular momentum. Until now, these layers required 
an ad-hoc treatment, such as performed by Brun et al. (1999) and Brun et al. (2002), who found that the mixing in the solar 
tachocline contributes indeed to the depletion of lithium, and that it modifies sufficiently the structure of the star to be detected 
through helioseismology. 

From a more technical point of view, we have written here the heat equation ( 110111021 . which determines the meridional 
circulation, as a relaxation equation for the temperature fluctuation ^F/. This breaks the apparent symmetry between the circulation 
driven by the differential rotation profile (hence by *F;) and the counteracting flow (Mestel's //-current) which is induced by the 
inhomogeneities in chemical composition (A/), since now the expression giving the circulation velocity involves only the first 
derivative in A;, while keeping the second derivative of 'F/. A benefit of that choice is that it renders this differential equation less 
sensitive, numerically, to the steep composition gradients which arise in the course of stellar evolution. 

One weakness of the modelization of rotational mixing remains the description of the turbulence which is generated by the 
differential rotation. Such turbulence is certainly anisotropic, due to the stable stratification, and it tends therefore to smooth the 
angular velocity and the chemical composition on horizontal surfaces. This problem is addressed in the companion paper (Mathis 
et al. 2004), where we discuss the prescriptions for the turbulent transport which are presently available. 

A major shortcoming of the present modelization of rotational mixing is that it predicts a rapidly spinning solar core, contrary 
to the findings of helioseismology. This proves that other physical processes contribute to the transport of angular momentum 
in solar-type stars, which are presently slow rotators because they have been spun down by their wind. Internal gravity waves, 
which are generated by turbulent convection at the interface with a convective region, are one possible cause. Their effect has 
been studied by Talon et al. (2002), and is now being introduced in stellar evolution codes (Talon & Charbonnel 2003). 

Magnetic stresses are another serious candidate for the angular momentum transport, as was pointed out already by Mestel 
(1953). The alternating dynamo field does not penetrate enough into the solar radiation zone to produce any impact on the 
rotation law, as was shown by Garaud (1999). But even a weak fossil field could enforce nearly uniform rotation, although the 
outcome depends sensitively on the poloidal field which is assumed, and on whether it penetrates into the differentially rotating 
convection zone, as illustrated by the calculations performed by Ruediger and Kitchanitov (1996, 1997) and by MacGregor 
and Charbonneau (1999). Therefore the problem must be treated consistently, taking into account the advection of the field by 
the meridional circulation, itself being modified by the field, and imposing proper boundary conditions. We shall address this 
problem in two forthcoming papers, one dealing with an axisymmetric field and the second with a non-axisymmetric field, as 
observed in oblique rotators. 
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Appendix A: Algebra related to the spherical harmonics 

The spherical harmonics are defined by: 

F,"'(6', if) = NfPY^ (cos 6) e'""f (A. 1) 

with the normahzation: 

2l+\il-\m\)\ 



NT-i-l)"^ 



^_^^^^^^^ and Br = — - 

V(2Z + l){2l - 1) V(2/ + 3)(2/+l) 



(A.2) 



47r (/-H|m|)! 
We recall that they obey the well-known differential equation: 

-de (sin OdeY",' {0, <f)) + -^dlYf (0, ^) = -/(/+ 1) F™ (6, <fi) . (A.3) 

sine ^ ' sin^e ^ 

In deriving certain equations in §5 we used the following recursion relations for m - 0: 



cos 0Y^{0) = A°Ff_j (0) + BlYf_^^{e) where AJ' = and B" = „^, „. J. ^7 ' (A-4) 



sin 0Y°i0) = C?dgY? , (0) - D°dgY^^A0) where C? = ^ — and D° = ^ . (A.5) 

' ' ' ' ' '^^ ' V(2Z-Hl)(2/-l) ' V(2/ + 3)(2/ + 1) 

cos0<9eFf(0) = E%yU0) + F%yUo) where £« = ^|J and Ff = [ (A.6) 

V(2/-i- 1)(2/- 1) -v(2( + 3)(2/ + 1) 

sin05eF,V) = G?}'° (0)-/7fF?,(0) where G? = —J^^i^^ and //? = — J^^l^^ . (A.7) 

' ' '^' ' ' ' ' V(2/ + 3)(2Z + 1) ' V(2/+l)(2/-l) 

The identities (IA.5> and ( IA.6> have been deduced from the two others JA.4HA.7l . with the help of JA.3I I. 

Appendix B: System of transport equations for 1-2,4 

In this appendix we give the explicit form of the transport equations for I - 2, A, with only the first term Q.2 being kept in the 
expansion of the rotation rate Q. 

B.1. Mean equations 
B.1.1. mean rotation rate 

p^ (r^Q) = ^dr (p/Qf/j) + 1^, (rVv..5,Q) (eq.[I3; (B.l) 

B.1.2. mean chemical composition 

P^-^ + ^dr [r^pTiUf''] = ^dr [MD„ + D,«)d,^] (eq.|36). (B.2) 
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B.2. System for 1=2 

B.2.1. meridional circulation 

We rewrite ( I101ll02l i in two first order equations as 

1 






^ad-^ 



S2 



where we have defined: 



S2 = -2 
1 






TiGp 

i + e, 



grav 



r^de 1-2 8- \ r 24- 2-2\ d f02 

g^dr\3 35 7 1I35 ' 3 / dr Uo 

^"•\''^<n '^^T L Di,\ 
r / X 1 M -id'i'j dlnll \ 



-QD.2 + - ( — Q2 - n|5,Q + —rQdr^2 

1 3 35 ^ / 35 



and: 



and where we have replaced the gravity fluctuation by j74> . /^ by ( I84t and /^^ by ( I85> . 
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(B.3) 



(B.4) 



(B.5) 



B.2.2. baroclinic relation 

From ( 146148 > we have: 



^A2 - (y^I'2 = = 



o 2 / 17 

7 3 35 7 35 



B.2.3. horizontal fluctuation of the molecular weight 
We apply (g^ to Z = 2: 

dK2 d\njl f/2 6 

-A2 = TT^f' - -tA,A2- 



f/f dt 



h/' r^- 



B.2.4. Poisson equation 
Likewise, for ( I60> : 



1 d^ 
r dr^ 



(r02) ■ 



;(;+!)_ 47rGdpo 



§0 dr 



47rG 
^0 



poa2 + :r (''Pobi) 
dr 



e.3. Sysfem for /=4 
B.3.1. horizontal shear 

d(r2Q2) 



dt 



2pQ.r 



■^^dr{pr^U2)-aU2 



-10pv/,n2 (eq-ED- 



B.3. 2. meridional circulation 
We proceed as for / = 2: 



f/4 






1 



:S4 



'ad ■ 



(B.6) 



(B.7) 



(B.8) 



(B.9) 



(B.IO) 
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where we have defined: 

,—2, 



-2 

1 
nGp 

1 



1 






f dr \35 



5H£(4no,)-i(|no.Ui|^ 



— (2QQ2 - '^ [^25,0 + Q5,Q2]) 



35 



r 20H: 

-OrJiA — 

3 3r 






+ e?rav r / \ n Af — /d*I'4 dlnu \ 



and: 



JU = //r5,-'I'4 - (1 - 5 + A-r) *P4 - (^ +Xi) A4 

B.3.3. baroclinic relation 

We restate J46l48t as 

(^A4 -5'¥a^ —^ \rQ.2drTi + rQ5,Q2 - 20021 . 
35 gL J 



B.3.4. horizontal fluctuation of the molecular weight 

dA4 din 77 f/4 20 ,_, 

— - ^A4 = -V, - -^D„A4 (eq.|33. 



B.3.5. Poisson equation 

i il f ^0 ) _ ^(^+1) 7 47rG dpo' 



47rG 



Pofl4 + :r ('■po^4) 

dr 



(eq.|60|l. 



(B.ll) 



(B.12) 



(B.13) 



(B.14) 



(B.15) 



r dr^ ^ '^'*> r^ ""* ^0 dr ^0 

These equations are ready to be implemented in a stellar evolution code, together with the boundary conditions discussed in 
§7. 



